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6.5 The Saddle Point Stokes Problem 



So far the matrix C has been diagonal — no trouble to invert. This section jumps to 
a fluid flow problem that is still linear (simpler than Navier- Stokes). But now 
represents the positive definite Laplacian (—A) in the continuous problem, and it is 
a finite difference or finite element discretization in the matrix problem. The matrix 
C^ 1 is sparse, but C itself involves a serious boundary value problem. 

In this case we don't reduce the computation all the way to A T CA\ Instead we 
stay with the sparse matrix whose blocks are C _1 , A, A T , and 0. This is an indefinite 
matrix, with that zero block on its diagonal. The solution v,p is a saddle point. 

This section also discusses the preconditioning of indefinite systems (not only 
Stokes). 



The Stokes Problem 

The unknowns in the two-dimensional Stokes problem are the velocity components 
in v = (vi(x,y),V2(x,y)) and the pressure p(x,y). The flow is incompressible so the 
velocity vector is constrained by divv = 0. Therefore A T is (minus) the divergence 
and A must be the gradient. Notice the different letters in our framework (u is 
changed to p, w is changed to v, and b is changed to /): 



pressure p(x, y) 



div v = 



A = gradient 



e = / — grad p 



v = Ce 



A T = —diver 



gence 



velocity v(x,y) 



Figure 6.15: The three-step A T CA framework for the Stokes problem. 



The real novelty is in e = C l v = (—Avi, — Av 2 )- The Stokes problem asks us to 
find a velocity vector v and a pressure p that solve —Av — f — gradp and divf = 0: 
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This describes slow viscous flow. It is not for aeronautics! The full Navier-Stokes 
equations have extra nonlinear terms from the motion of the underlying fluid. The 
Stokes problem could come into its own in biological applications, not for large-scale 
blood flow in the heart but for small-scale movements in capillaries (or even in cells). 
This linear problem is simpler than Navier-Stokes, and needs to be understood first. 
We keep it simple by omitting discussion of the boundary conditions. 

In the language of optimization, the pressure p(x, y) is the Lagrange multiplier 
that imposes the incompressibility constraint div v — when the energy is minimized: 



Minimize 



gradt>i| 2 + | gradt> 2 | 2 — 2f ■ gradf) dxdy subject to divt> = 



The Lagrangian L(v,p) adds J J 5L/5v = and 5L/5p = are exactly the Stokes 
equations (1). Chapter 6 will develop this primal-dual block form to find the saddle 
point of L — here we only need to know A and C. 

The key point is that we do not eliminate v to reach K = A T CA, because C is 
the inverse of the Laplacian. It is diagonal only in frequency space, not in physical 
x, y space. The Stokes problem stays in its block form (1) with two unknowns v and 
p: 

Block form (continuous) 



'C- 1 A' 
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—A grad 
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In the discrete case, that block matrix will be symmetric but indefinite (as usual). 
It will have positive and also negative eigenvalues and pivots (positive pivots from 
C _1 , negative pivots from — A T CA). Finite elements for velocity v and pressure p 
have to be carefully matched to ensure that the discrete problem has a good solution. 
The approximating polynomials are often one degree higher for v than for p. 



The Inf-Sup Condition 

Restricting a positive definite operator if to a subspace (by a projection P) never 
destroys positive definiteness. All eigenvalues of PKP T lie safely between X m i n (K) 
and A max (if). But an indefinite operator has a negative A m i n , so this interval includes 
zero: not safe! We must take particular care that the block matrix is safely invertible. 

These mixed problems, or saddle point problems, are partly positive and 
partly negative. The new requirement for a bounded inverse is an inf-sup condition 
(also known as a Babuska-Brezzi condition): For every p there must be a v so 
that 

v T Ap > f3 Vv T C- 1 v ^fp?p~ for a fixed (3 > . (3) 
This "compatibility" of velocities and pressures produces ||(A T CA) _1 || < 

Condition (3) will immediately fail if there is a nonzero pressure p with Ap = 0. 
In the discrete case, this means that A must have full column rank. Otherwise the 
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last columns of the block matrix in (2) are not independent, and the block matrix is 
not invertible. In particular, the rectangular A cannot be short and wide. 

More exactly, the lower bound (3) on v T Ap requires that a nonzero pressure p 
may not be orthogonal to every A T v. The space of all A T v 's must have at least 
the dimension of the space of p's. Thinking of A T as divergence, this indicates why 
the velocity finite elements are often one polynomial degree above the pressure finite 
elements. But dimensions of spaces are not sufficient to confirm that (3) is actually 
true. Each choice of finite elements for v and p requires its own analysis. 

Now we explain the highlighted statement after equation (3): A T CA is bounded 
below by f3 2 and ||(.A T CL4) _1 || < 1//9 2 when the inf-sup condition holds. Begin by 
introducing w = C~ 1//2 t> . The continuous and discrete Laplacians, both called C -1 , 
have positive definite square roots C~ l l 2 (never to be computed). Then (3) reads: 



^ w T CV 2 Ap 

box every p there must be a if so that — 

\\w\\ 



> PWpW • (4) 



This is where we take the maximum over all w (the sup). On the left, the maximum 
of w T z/||w|| is exactly ||z||. This is attained when w = z (then the cosine is one). 
Since (4) has z = C 1 ^ 2 Ap, this best choice of w reduces (3) to 

For every p, \\z\\ = \\C 1/2 Ap\\ > (3\\p\\ . (5) 

Square both sides, and minimize their ratio over all p (the inf). The inequality 

\\C 1/2 Ap\\ 2 > p 2 \\p\\ 2 or p T A T CAp>(3 2 p T p (6) 

is tightest at the lowest eigenvector and eigenvalue, where A T CAp = j3 2 p. The largest 
eigenvalue of (A^-CA)" 1 is 1//3 2 . The matrix is symmetric so its norm is also 1//3 2 . 

A key goal in replacing || (A^CA) -1 1| < 1/ (3 2 by the inf-sup condition is to work 
with C~ l and not C. We can apply (3) to the continuous problem (it holds, so Stokes 
is well-posed). We can also test (3) on subspaces Vh and Ph of finite element trial 
functions. It is a compatibility condition: for every pressure pn in Ph, there must be 
a velocity Vh in Vh satisfying (3). If and when this is established, with a bound flh 
that stays away from zero, the mixed finite element method will be stable. 

Correction to p T p: Finite elements are functions ph(x,y), not just vectors v and 
p. So inner products and norms of Vh and ph come from integrals and not sums. The 
L 2 norm of a pressure trial function ph = YlPj^jfaiV) ^ s connected to the discrete 
coefficient vector p through a positive definite "pressure mass matrix" Q: 

\\Ph\\ 2 = J J iphf dx dy = ^2^PiPj J j Mj dx dy = p T Qp . (7) 

The correct inf-sup condition (3) changes p T p to p T Qp. When satisfied for the finite 
elements 3, 5, 8 in the list below, the inequality (6) involves p T Qp: 

Correction to (6) p T A T CAp > (3 2 p T Qp (8) 

This makes Q (or even its diagonal part) a simple and useful candidate as precondi- 
tioner. Not the best, but it is directly available and it scales the problem correctly. 
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Testing the Inf-Sup Condition 

The inf-sup condition can be delicate. The books by Brezzi-Fortin and Elman-Silvester- 
Wathen are excellent references. From the latter book we report on typical results, 
writing Pq,Pi,P 2 for constant, linear, and quadratic polynomials on triangles, and 
Q01Q11Q2 for constant, bilinear, and biquadratic elements on quadrilaterals: 



1. Velocities in Pi, pressures in P : failure 

2. Velocities in P 1; pressures in P x : failure 

3. Velocities in P 2 , pressures in P\. success 

4. Velocities in Qi, pressures in Q\\ failure 

5. Velocities in Q 2 , pressures in Qi. success (Figure 6.17) 





Figure 6.16: Q 2 velocities • and Qi pressures : failure on one square, success on 
two. 



For Q 2 — Qi on one square, A is 2 by 4 (four pressures, two equations div v = 
at the center). There will surely be non-constant solutions to Ap = 0: failure. 

For two squares, A is 6 by 6 (six pressures, divv = at three internal nodes). 
Only a constant pressure solves Ap = 0, and this hydrostatic solution is removed 
by the boundary conditions (like grounding a node, here a better choice would be 
J J pdx dy = 0). A related "serendipity" element is also successful, by removing the 
node inside the square and the x 2 y 2 term in the velocities. That leaves 8 nodes 
matching 8 terms in the polynomial. 

The list of element pairs also includes pressures P_i and Q_i that are not contin- 
uous between elements (this is permitted since delta functions don't appear): 



6. Velocities in P 2 , pressures in P_ x : failure 

7. Velocities in Q 2 , pressures in Q-\. failure 

8. Velocities in Q 2 , pressures in P_l: success 

9. Velocities in Qi, pressures in Q : failure (Figure 6.18) 
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In each square 
v — a + bx + cy + dxy in Qi 
p = constant in Q 



Figure 6.17: Bilinear velocities •, constant pressures O (unstable checkerboard 
mode). 

The failure of Qi — Qo is especially important. Qo means a constant pressure in each 
square. This is the simplest pair of elements that are conforming (no delta functions 
in gradv because Q\ elements a + bx + cy + dxy are continuous across edges). 

Four squares are not stable because that pressure vector p — (1, —1, 1, —1) satisfies 
Ap = 0. More squares don't help — this checkerboard pattern with "white squares 
= 1" and "black squares = —1" continues to be a spurious mode satisfying Ap = 0. 
The columns of A, even after p = constant is removed, are still not independent. 

This Q\— Qo element is so simple and potentially useful that it is often stabilized by 
relaxing the incompressibility equation A T v = 0. We insert a new anti-checkerboard 
matrix — aE, exactly chosen to eliminate (1, —1, 1, —1) from the nullspace! 



Stabilized The block matrix becomes 



c- 1 

A T 



A 

-aE 



with E 



1 -1 

•1 1 

1 -1 

•1 1 



1 -1 

•1 1 

1 -1 

•1 1 



Our zero block is replaced by the negative semidefinite matrix —aE. Elimination 
leaves —A T CA — aE in that block, which is more strongly negative. The Qi — Qo 
construction is now uniformly invertible, since E times the checkerboard p is not zero. 

In this special case we know the guilty p. In other cases E is constructed with only 
the constant pressure in its nullspace (which isn't harmful). The new aE is taken 
small enough so that its positive eigenvalues roughly match those of A T CA. The zero 
eigenvalue of A T CA (produced by the checkerboard pressure) becomes positive. 



Solving the Discrete Saddle Point Problem 

Don't forget that this is the fundamental problem of scientific computing! 

or A T CAu = f - A T Cb. (9) 



~c- 1 


A' 




w 




'b~ 


A T 







u 







In this section, w has become v and u is p. In earlier sections, when C was a 
diagonal matrix (small-block diagonal, in a multidimensional problem), the reduction 
to A T CA was very reasonable. Now the Stokes problem illustrates that C~ l can be 
sparse (coming from finite differences or finite elements) while C is full. In that case 
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we may prefer to compute with the sparse block form, even though its matrix is not 
at all positive definite like A T CA. 

A direct method (elimination) is appropriate up to n = 10 4 or 10 5 . It will be used 
in most problems! Of course an indefinite matrix may require row exchanges to avoid 
small pivots. A good code will decide that without our intervention. 

The pivoting can be organized using 2 by 2 blocks, or we may reorder the unknowns 

as in Section to keep the L and U factors as sparse as possible. If we eliminate 

in the original order, subtracting A T C times one block row from the other, then 
— K = —A T CA will appear as the Schur complement in the lower right block. 

For very large systems, a direct method demands too much computation (and 
memory). We'll look briefly at iterative methods, remembering that the conjugate 
gradient method expects positive definiteness. It could be used with A T CA, but 
it is not safe for the indefinite problem. Each CG step applies A T CA (possibly 
preconditioned) to a vector. The matrix- vector multiplications Ap and A T v are fast. 
In between, multiplying by C amounts to solving the system C~ x v = Ap. In practice 
this will be two linear systems for the velocity components v\ and t>2, with the same 
Laplace matrix — two "Poisson solves" at each step. This is not impossible but some 
inner iteration like multigrid is likely to be used. 

A number of iterative methods compute both v and p, based on the block form: 

1. The Uzawa method has C^v k+1 = f — Ap k and p k+1 = p k + aA T v k+1 

2. Penalty methods are described in Chapter on optimization 

3. The augmented Lagrangian method is presented in detail by Glowinski-Le 
Tallec. 

These and more are discussed by Quarteroni-Valli for the Stokes problem. They all 
involve the "Poisson solve" associated with C. We turn instead to a closer look at 
preconditioners for the general saddle point problem, whether it arises in fluid flow or 
optimization or network analysis or elsewhere. 

Preconditioning Saddle Point Problems 

The basic methods for large sparse indefinite problems are MINRES and GMRES 
(symmetric and nonsymmetric). Our focus here is on constructing a good precondi- 
tioner. 

When we only know the matrix entries in C~ l and A, the code will have to 
construct a good preconditioner "blindly" from that information. Recall the funda- 
mental dividing line mentioned in the preface, to be given a matrix or given a problem. 
When we know the underlying problem (like Stokes), we can devise a preconditioner 
by hand. There is a similar separation between algebraic multigrid and geometric 
multigrid (the computer thinks or we do). 
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Look first at the saddle point matrix M and its block factorization M = LDL 



T, 



M 
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(10) 



'd- 1 a' 




c- 1 


P 3 = 


d~ l a 


_A T 0_ 


P 2 = 


K 


K 



M~ l clearly involves C and (A T CA)~ 1 from the inverse of that middle matrix. Our 
preconditioners will keep this m, n block structure. To come close to M _1 , as precon- 
ditioned are intended to do, we could approximate the blocks of M and invert or we 
could directly approximate C and (A T CA)~ l . 

Three approximations to M are suggested by Elman-Silvester-Wathen: 



Preconditioners Pi 



The matrix A (often relatively simple) is kept unchanged. C approximates the Poisson 
solver C, and K approximates K = A T CA. We must choose C and K. 

A key point about P2 and P3: If C = C and K = K } then P 2 ~ 1 M has only three 
different eigenvalues and P^ l M has only the eigenvalue A = 1. In this case. In that 
case MINRES and GMRES converge in three or two steps. So those changes from Pi 
add very few iterations. Wathen begins with simple suggestions for the preconditioner 

1. Replace the Poisson solver C by a single multigrid cycle C 

2. Replace K~ l by four conjugate gradient steps preconditioned by Q or diag(Q). 



The multigrid cycle applies a basic iteration on fine and also coarse meshes, to reduce 

both high and low frequencies in the error residual (Section ). The matrix Q 

that preconditions K is the pressure mass matrix in (7). Its energy p T Qp gives a 
lower (and also upper) bound to p T Kp. Usually diag(Q) also has this "equispectral 
property", scaling the system correctly even if the approximation to K is not really 
close. 

C and K are defined algorithmically by multigrid and preconditioned conjugate 
gradient cycles. This is simpler than trying to express them by a matrix. 



Nonsymmetric Problems and Model Reduction 

We move now toward problems more general than Stokes. They might be associated 
with the Navier-Stokes equations (div v is still zero) and its linearizations. In that case 
C will be unsymmetric, from the first derivative terms in the linearization. Convection 
joins diffusion (also in many problems from chemical engineering). 

A key point is the relative strength of those two terms, inertia versus viscosity. It is 
measured by the Reynolds number Re = (density) (velocity) (distance) /(viscosity). 
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Stokes flow is the limiting case Re = for high viscosity and low velocity. The high 
speed Euler equation dv/dt = — gradp + F is at the other extreme Re — » oo. 

Beyond computational fluid dynamics are other differential equations with their 
own constraints, and then all the algebraic problems of constrained optimization. 
Our discussion is moving from "given a problem" toward "given a matrix." The 
constraints make it a saddle point matrix — possibly unsymmetric, with A T changed 
to B. We look to preconditioners of the same general forms, as a guide. 

Always the key is an approximate inverse to the Schur complement —K = —BCA: 



Elimination reduces 



'C- 1 A' 




to 


B 



c- 1 





A 
-BCA 



c- 1 





A 

-K 



Normally K will be dense. If it is also large, we look at iterative methods. Multigrid 
for C will become algebraic multigrid [Briggs-McCormick-van Hensen]. The approx- 
imation is constructed directly from the matrix, not using the differential equation 
that may be in the background. 

For K = BCA, a "black box" approach to preconditioning is suggested in [ ]: 

Reduce C to a smaller model matrix C, satisfying 



AC ~ CA and then K — BAC « BCA. 



:i2) 



Notice that C is m by m (large) while C is n by n. The model reduction CA = AC 
will not be exactly solvable (more equations than unknowns). The meaning of ~ is to 
be decided by the user! We may apply weighted least squares to compute C~ l from 
AC~ l ~ C~ 1 A, since the saddle point problem gives us C~ l and not C. 

I expect new ideas will keep coming, for this fundamental indefinite problem. 



Ad vect ion-Diffusion Problems 

The Stokes problem involves the Laplacian — Av. The corresponding matrix C~ l is 
"essentially symmetric." Perfect symmetry holds in the interior and could only be 
lost at the boundary (it shouldn't be). Now, linearizing the v ■ gradf term in the 
Navier-Stokes equation produces a first derivative "advection" term, and symmetry 
is definitely lost. 



